#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <ctype.h>
#include <float.h>
#include <string.h>
#include <stdarg.h>
#include <limits.h>
#include <locale.h>
#include "svm.h"
int libsvm_version = LIBSVM_VERSION;
typedef float Qfloat;
typedef signed char schar;
#ifndef min
template <class T> static inline T min(T x,T y) { return (x<y)?x:y; }
#endif
#ifndef max
template <class T> static inline T max(T x,T y) { return (x>y)?x:y; }
#endif
template <class T> static inline void swap(T& x, T& y) { T t=x; x=y; y=t; }
template <class S, class T> static inline void clone(T*& dst, S* src, int n)
{
  dst = new T[n];
  memcpy((void *)dst,(void *)src,sizeof(T)*n);
}
static inline double powi(double base, int times)
{
  double tmp = base, ret = 1.0;

  for(int t=times; t>0; t/=2)
    {
      if(t%2==1) ret*=tmp;
      tmp = tmp * tmp;
    }
  return ret;
}
#define INF HUGE_VAL
#define TAU 1e-12
#define Malloc(type,n) (type *)malloc((n)*sizeof(type))

static void print_string_stdout(const char *s)
{
  fputs(s,stdout);
  fflush(stdout);
}
static void (*svm_print_string) (const char *) = &print_string_stdout;
#if 1
static void info(const char *fmt,...)
{
  char buf[BUFSIZ];
  va_list ap;
  va_start(ap,fmt);
  vsprintf(buf,fmt,ap);
  va_end(ap);
  (*svm_print_string)(buf);
}
#else
static void info(const char *fmt,...) {}
#endif

//
// Kernel Cache
//
// l is the number of total data items
// size is the cache size limit in bytes
//
class Cache
{
public:
  Cache(int l,long int size);
  ~Cache();

  // request data [0,len)
  // return some position p where [p,len) need to be filled
  // (p >= len if nothing needs to be filled)
  int get_data(const int index, Qfloat **data, int len);
  void swap_index(int i, int j);	
private:
  int l; // 样本总数。
  long int size; // 指定的全部内存，据说用Mb做单位。
  struct head_t
  {
    head_t *prev, *next;	// a circular list
    Qfloat *data;
    int len;		// data[0,len) is cached in this entry
  };

  head_t *head; // 记录申请来的内存，双向链表记录。
  head_t lru_head; // 双向链表头
  void lru_delete(head_t *h);
  void lru_insert(head_t *h);
};
// 根据样本数L申请L个head_t的空间。
Cache::Cache(int l_,long int size_):l(l_),size(size_)
{
  head = (head_t *)calloc(l,sizeof(head_t));	// initialized to 0 
  size /= sizeof(Qfloat); // 单位化内存记录单位 
  size -= l * sizeof(head_t) / sizeof(Qfloat); // 扣除已申请到的内存大小
  size = max(size, 2 * (long int) l);	// cache must be large enough for two columns
  lru_head.next = lru_head.prev = &lru_head;
}

Cache::~Cache()
{
  for(head_t *h = lru_head.next; h != &lru_head; h=h->next)
    free(h->data);
  free(head);
}

void Cache::lru_delete(head_t *h)
{
  // delete from current location
  h->prev->next = h->next;
  h->next->prev = h->prev;
}

void Cache::lru_insert(head_t *h)
{
  // insert to last position
  h->next = &lru_head;
  h->prev = lru_head.prev;
  h->prev->next = h;
  h->next->prev = h;
}
// 该函数保证head[index]中至少有len个float的内存，
int Cache::get_data(const int index, Qfloat **data, int len)
{
  head_t *h = &head[index];
  if(h->len) lru_delete(h);
  int more = len - h->len;

  if(more > 0)
    {
      // free old space
      while(size < more)
	{
	  head_t *old = lru_head.next;
	  lru_delete(old);
	  free(old->data);
	  size += old->len;
	  old->data = 0;
	  old->len = 0;
	}

      // allocate new space
      h->data = (Qfloat *)realloc(h->data,sizeof(Qfloat)*len);
      size -= more;
      swap(h->len,len);
    }

  lru_insert(h);
  *data = h->data;
  return len;
}

void Cache::swap_index(int i, int j)
{
  if(i==j) return;

  if(head[i].len) lru_delete(&head[i]);
  if(head[j].len) lru_delete(&head[j]);
  swap(head[i].data,head[j].data);
  swap(head[i].len,head[j].len);
  if(head[i].len) lru_insert(&head[i]);
  if(head[j].len) lru_insert(&head[j]);

  if(i>j) swap(i,j);
  for(head_t *h = lru_head.next; h!=&lru_head; h=h->next)
    {
      if(h->len > i)
	{
	  if(h->len > j)
	    swap(h->data[i],h->data[j]);
	  else
	    {
	      // give up
	      lru_delete(h);
	      free(h->data);
	      size += h->len;
	      h->data = 0;
	      h->len = 0;
	    }
	}
    }
}

//
// Kernel evaluation
//
// the static method k_function is for doing single kernel evaluation
// the constructor of Kernel prepares to calculate the l*l kernel matrix
// the member function get_Q is for getting one column from the Q Matrix
//
class QMatrix {
public:
  virtual Qfloat *get_Q(int column, int len) const = 0; // 返回第column列向量
  virtual double *get_QD() const = 0;
  virtual void swap_index(int i, int j) const = 0;
  virtual ~QMatrix() {}
};

class Kernel: public QMatrix {
public:
  Kernel(int l, svm_node * const * x, const svm_parameter& param);
  virtual ~Kernel();

  static double k_function(const svm_node *x, const svm_node *y,
			   const svm_parameter& param);
  virtual Qfloat *get_Q(int column, int len) const = 0; // 纯虚函数
  virtual double *get_QD() const = 0;
  virtual void swap_index(int i, int j) const	// no so const...
  {
    swap(x[i],x[j]);
    if(x_square) swap(x_square[i],x_square[j]);
  }
protected:

  double (Kernel::*kernel_function)(int i, int j) const;

private:
  const svm_node **x; // 指向样本数据，每次数据传入时通过克隆函数来实现，完全重新分配内存。
  double *x_square; // 使用RBF核才使用

  // svm_parameter
  const int kernel_type; 
  const int degree;
  const double gamma;
  const double coef0;

  static double dot(const svm_node *px, const svm_node *py); // 点乘
  double kernel_linear(int i, int j) const 
  {
    return dot(x[i],x[j]);
  }
  double kernel_poly(int i, int j) const
  {
    return powi(gamma*dot(x[i],x[j])+coef0,degree);
  }
  double kernel_rbf(int i, int j) const
  {
    return exp(-gamma*(x_square[i]+x_square[j]-2*dot(x[i],x[j])));
  }
  // tanh(x) = (e^x-e^(-x)) /(e^x+e^(-x)) ;
  // 计算双曲正切函数
  double kernel_sigmoid(int i, int j) const
  {
    return tanh(gamma*dot(x[i],x[j])+coef0);
  }
  double kernel_precomputed(int i, int j) const
  {
    return x[i][(int)(x[j][0].value)].value;
  }
};

Kernel::Kernel(int l, svm_node * const * x_, const svm_parameter& param)
  :kernel_type(param.kernel_type), degree(param.degree),
   gamma(param.gamma), coef0(param.coef0)
{
  // 根据核函数类型，初始化核函数指针
  switch(kernel_type)
    {
    case LINEAR:
      kernel_function = &Kernel::kernel_linear;
      break;
    case POLY:
      kernel_function = &Kernel::kernel_poly;
      break;
    case RBF:
      kernel_function = &Kernel::kernel_rbf;
      break;
    case SIGMOID:
      kernel_function = &Kernel::kernel_sigmoid;
      break;
    case PRECOMPUTED:
      kernel_function = &Kernel::kernel_precomputed;
      break;
    }

  clone(x,x_,l);

  if(kernel_type == RBF)
    {
      x_square = new double[l];
      for(int i=0;i<l;i++)
	x_square[i] = dot(x[i],x[i]);
    }
  else
    x_square = 0;
}

Kernel::~Kernel()
{
  delete[] x;
  delete[] x_square;
}

double Kernel::dot(const svm_node *px, const svm_node *py)
{
  double sum = 0;
  while(px->index != -1 && py->index != -1)
    {
      if(px->index == py->index)
	{
	  sum += px->value * py->value;
	  ++px;
	  ++py;
	}
      else
	{
	  if(px->index > py->index)
	    ++py;
	  else
	    ++px;
	}			
    }
  return sum;
}
// 核函数，只有在预报的时候才用到
double Kernel::k_function(const svm_node *x, const svm_node *y,
			  const svm_parameter& param)
{
  switch(param.kernel_type)
    {
    case LINEAR:
      return dot(x,y);
    case POLY:
      return powi(param.gamma*dot(x,y)+param.coef0,param.degree);
    case RBF:
      {
	double sum = 0;
	while(x->index != -1 && y->index !=-1)
	  {
	    if(x->index == y->index)
	      {
		double d = x->value - y->value;
		sum += d*d;
		++x;
		++y;
	      }
	    else
	      {
		if(x->index > y->index)
		  {	
		    sum += y->value * y->value;
		    ++y;
		  }
		else
		  {
		    sum += x->value * x->value;
		    ++x;
		  }
	      }
	  }

	while(x->index != -1)
	  {
	    sum += x->value * x->value;
	    ++x;
	  }

	while(y->index != -1)
	  {
	    sum += y->value * y->value;
	    ++y;
	  }
			
	return exp(-param.gamma*sum);
      }
    case SIGMOID:
      return tanh(param.gamma*dot(x,y)+param.coef0);
    case PRECOMPUTED:  //x: test (validation), y: SV
      return x[(int)(y->value)].value;
    default:
      return 0;  // Unreachable 
    }
}

// An SMO algorithm in Fan et al., JMLR 6(2005), p. 1889--1918
// Solves:
//
//	min 0.5(\alpha^T Q \alpha) + p^T \alpha
//              p is the vector of all ones
//		y^T \alpha = \delta
//		y_i = +1 or -1
//		0 <= alpha_i <= Cp for y_i = 1
//		0 <= alpha_i <= Cn for y_i = -1
//              Q is an l by l symmetric matrix with Qij = Yi*Yj*K(Xi,Xj) and K(Xi,Xj) is kernel function.
// Given:
//
//	Q, p, y, Cp, Cn, and an initial feasible point \alpha
//	l is the size of vectors and matrices
//	eps is the stopping tolerance
//
// solution will be put in \alpha, objective value will be put in obj
//
class Solver {
public:
  Solver() {};
  virtual ~Solver() {};

  struct SolutionInfo {
    double obj;
    double rho;
    double upper_bound_p;
    double upper_bound_n;
    double r;	// for Solver_NU
  };

  void Solve(int l, const QMatrix& Q, const double *p_, const schar *y_,
	     double *alpha_, double Cp, double Cn, double eps,
	     SolutionInfo* si, int shrinking);
protected:
  int active_size; // 计算时实际参加运算的样本数目，经过shrink处理后，该数目会小于全部样本总数
  schar *y; // 样本所属类别，该值为+1/-1 
  double *G;		// gradient of objective function 存储的是dilta用来统计b的上下界 找到ktt的反例 进行迭代的。
  // Q*alpha[i]+p[i] = G[i] 
  enum { LOWER_BOUND, UPPER_BOUND, FREE }; // a[i] = 0 , YiUi>1的点，正确分类，且不在分类线上的点
  // a[i] = C , YiUi<1的点，在两条分类线间的点 
  // 0 =< a[i] <= C YiUi = 1 的点，在分类线上的点
  char *alpha_status;	// LOWER_BOUND , UPPER_BOUND, FREE 
  double *alpha; // a[i] 
  const QMatrix *Q; // 指定核。核函数和Solver相互结合，可以产生多种SVC，SVR 
  const double *QD; // Q对角矩阵的值
  double eps; // 误差
  double Cp,Cn; // a[i]的上限 惩罚因子
  double *p; // 用来计算gradient
  int *active_set;
  double *G_bar;		// gradient, if we treat free variables as 0
  // G_bar[i] = C*sum(Q[k][i]) k = 1,...,l when a[i] == C
  // 此向量用在做shrink的时候减少梯度计算亮 
  // G[i] = G_bar[i] + sum(Q[k][i]*a[i]) when a[i] != C 
  int l; // 样本总数
  bool unshrink;	// XXX
  // 返回对应于样本的c。设置不同的Cp和Cn是为了处理数据的不平衡 有时Cp=Cn
  double get_C(int i) 
  {
    return (y[i] > 0)? Cp : Cn;
  }
  // 做边界处理
  void update_alpha_status(int i)
  {
    if(alpha[i] >= get_C(i))
      alpha_status[i] = UPPER_BOUND;
    else if(alpha[i] <= 0)
      alpha_status[i] = LOWER_BOUND;
    else alpha_status[i] = FREE;
  }
  bool is_upper_bound(int i) { return alpha_status[i] == UPPER_BOUND; }
  bool is_lower_bound(int i) { return alpha_status[i] == LOWER_BOUND; }
  bool is_free(int i) { return alpha_status[i] == FREE; }
  void swap_index(int i, int j);
  void reconstruct_gradient();
  virtual int select_working_set(int &i, int &j);
  virtual double calculate_rho();
  virtual void do_shrinking();
private:
  bool be_shrunk(int i, double Gmax1, double Gmax2);	
};
// 完全交换样本i和j的内容，包括所申请的内存的地址
void Solver::swap_index(int i, int j)
{
  Q->swap_index(i,j); // 交换核矩阵
  swap(y[i],y[j]); // 交换标签
  swap(G[i],G[j]); // 交换梯度 
  swap(alpha_status[i],alpha_status[j]); // 交换alpha状态
  swap(alpha[i],alpha[j]); // 交换alpha的值
  swap(p[i],p[j]);
  swap(active_set[i],active_set[j]);
  swap(G_bar[i],G_bar[j]);
}
/* 
   重新计算梯度。
   
*/
void Solver::reconstruct_gradient()
{
  // reconstruct inactive elements of G from G_bar and free variables

  if(active_size == l) return;

  int i,j;
  int nr_free = 0;

  for(j=active_size;j<l;j++)
    G[j] = G_bar[j] + p[j];

  for(j=0;j<active_size;j++)
    if(is_free(j))
      nr_free++;

  if(2*nr_free < active_size)
    info("\nWARNING: using -h 0 may be faster\n");

  if (nr_free*l > 2*active_size*(l-active_size))
    {
      for(i=active_size;i<l;i++)
	{
	  const Qfloat *Q_i = Q->get_Q(i,active_size);
	  for(j=0;j<active_size;j++)
	    if(is_free(j))
	      G[i] += alpha[j] * Q_i[j];
	}
    }
  else
    {
      for(i=0;i<active_size;i++)
	if(is_free(i))
	  {
	    const Qfloat *Q_i = Q->get_Q(i,l);
	    double alpha_i = alpha[i];
	    for(j=active_size;j<l;j++)
	      G[j] += alpha_i * Q_i[j];
	  }
    }
}

void Solver::Solve(int l, const QMatrix& Q, const double *p_, const schar *y_,
		   double *alpha_, double Cp, double Cn, double eps,
		   SolutionInfo* si, int shrinking)
{
  this->l = l;
  this->Q = &Q;
  QD=Q.get_QD();
  clone(p, p_,l);
  clone(y, y_,l);
  clone(alpha,alpha_,l);
  this->Cp = Cp;
  this->Cn = Cn;
  this->eps = eps;
  unshrink = false;

  // initialize alpha_status
  {
    alpha_status = new char[l];
    for(int i=0;i<l;i++)
      update_alpha_status(i);
  }

  // initialize active set (for shrinking)
  {
    active_set = new int[l];
    for(int i=0;i<l;i++)
      active_set[i] = i;
    active_size = l;
  }

  // initialize gradient
  {
    G = new double[l];
    G_bar = new double[l];
    int i;
    for(i=0;i<l;i++)
      {
	G[i] = p[i];
	G_bar[i] = 0;
      }
    for(i=0;i<l;i++)
      if(!is_lower_bound(i)) // 计算 alpha > 0 的 梯度函数
	{
	  const Qfloat *Q_i = Q.get_Q(i,l); // 得到第i列的Q 
	  double alpha_i = alpha[i]; // 得到第i个向量的alpha
	  int j;
	  for(j=0;j<l;j++)
	    G[j] += alpha_i*Q_i[j];  // 计算出梯度 G = Qa+e ; G[j] = Sum( Q[j][i]*a[i])when(0<=i<l) + 1 ; 
	  if(is_upper_bound(i)) 
	    for(j=0;j<l;j++)
	      G_bar[j] += get_C(i) * Q_i[j]; // G_bar 保存的是在a=c的时候的C*Q;
	}
  }

  // optimization step

  int iter = 0;
  int max_iter = max(10000000, l>INT_MAX/100 ? INT_MAX : 100*l);
  int counter = min(l,1000)+1;
	
  while(iter < max_iter)
    {
      // show progress and do shrinking

      if(--counter == 0)
	{
	  counter = min(l,1000);
	  if(shrinking) do_shrinking();
	  info(".");
	}

      int i,j;
      if(select_working_set(i,j)!=0)
	{
	  // reconstruct the whole gradient
	  reconstruct_gradient();
	  // reset active set size and check
	  active_size = l;
	  info("*");
	  if(select_working_set(i,j)!=0)
	    break;
	  else
	    counter = 1;	// do shrinking next iteration
	}
		
      ++iter;

      // update alpha[i] and alpha[j], handle bounds carefully
		
      const Qfloat *Q_i = Q.get_Q(i,active_size);
      const Qfloat *Q_j = Q.get_Q(j,active_size);

      double C_i = get_C(i);
      double C_j = get_C(j);

      double old_alpha_i = alpha[i];
      double old_alpha_j = alpha[j];

      if(y[i]!=y[j])
	{
	  double quad_coef = QD[i]+QD[j]+2*Q_i[j]; // Kii+Kjj-2*Kji Q为 Yi*Yj*Kij 当Yi!=Yj Qji = -Kji ;
	  if (quad_coef <= 0)
	    quad_coef = TAU; // 如果quad_coef <= 0 那么极值在边界处
	  
	  // 计算出更新的值 alpha_i_new = alpha_i_old + Yj(Ei-Ej);
	  // Ei = Sum(alpha[i]*Yj*K[j][i])when(0=<j<l) - Yi ;
	  // 通过把Yj与E结合，delta等于上述的等式
	  double delta = (-G[i]-G[j])/quad_coef;  
	  
	  // Yi*alpha[i]+Yj*alpha[j] = constant 
	  // 故在Yi Yj异号时候 alpha[i]-alpha[j] = constant ; 
	  double diff = alpha[i] - alpha[j]; 
	  // alpha[i]*Yi + alpha[j]*Yj = alpha_old[i]*Yi + alpha_old[j]*Yj = constant ; alpha[i],alpha[h] is variable and the other alpha is constant ; so alpha[i]*Yi = alpha_old[i]*Yi + Yj * (alpha_old[j]-alpha[j]) ;
	  // alpha[j] = alpha_old[j] + delta ;
	  // alpha[i] = alpha_old[i] + Yi*Yj*(alpha_old[j]-alpha[j]) ;
	  //          = alpha_old[i] + delta ;
	  alpha[i] += delta;
	  alpha[j] += delta;
		

	  // 对结果进行截断
	  /*
	    alpha = H if alpha >= H ;
	            alpha if alpha < H and alpha > L ;
		    L if alpha <= L ;
	    L = max(0,ai-aj) H = min(C,C+ai-aj) ; when Yi != Yj ;
	    L = max(0,ai-aj-C) H = min(C,ai+ajf ; when Yi == Yj ;
	   */
	  if(diff > 0) // diff = ai - aj  L = diff H = C ;
	    {
	      if(alpha[j] < 0) // aj < 0 when ai <= L ;
		{
		  alpha[j] = 0; //  aj = aj - (ai - ai + aj) ;
		  alpha[i] = diff; // ai = L = ai - aj ;
		}
	    }
	  else // L = 0 ; H = C + ai-aj;
	    {
	      if(alpha[i] < 0) 
		{
		  alpha[i] = 0;
		  alpha[j] = -diff; // aj = aj - (ai-0) = aj -ai = -diff ;
		}
	    }
	  if(diff > C_i - C_j) //  
	    {
	      if(alpha[i] > C_i) // alpha > H 
		{
		  alpha[i] = C_i; 
		  alpha[j] = C_i - diff;
		}
	    }
	  else
	    {
	      if(alpha[j] > C_j)
		{
		  alpha[j] = C_j; 
		  alpha[i] = C_j + diff;
		}
	    }
	}
      else
	{
	  double quad_coef = QD[i]+QD[j]-2*Q_i[j];
	  if (quad_coef <= 0)
	    quad_coef = TAU;
	  double delta = (G[i]-G[j])/quad_coef;
	  double sum = alpha[i] + alpha[j];
	  alpha[i] -= delta;
	  alpha[j] += delta;

	  if(sum > C_i)
	    {
	      if(alpha[i] > C_i)
		{
		  alpha[i] = C_i;
		  alpha[j] = sum - C_i;
		}
	    }
	  else
	    {
	      if(alpha[j] < 0)
		{
		  alpha[j] = 0;
		  alpha[i] = sum;
		}
	    }
	  if(sum > C_j)
	    {
	      if(alpha[j] > C_j)
		{
		  alpha[j] = C_j;
		  alpha[i] = sum - C_j;
		}
	    }
	  else
	    {
	      if(alpha[i] < 0)
		{
		  alpha[i] = 0;
		  alpha[j] = sum;
		}
	    }
	}

      // update G

      double delta_alpha_i = alpha[i] - old_alpha_i;
      double delta_alpha_j = alpha[j] - old_alpha_j;
      
      for(int k=0;k<active_size;k++)
	{
	  G[k] += Q_i[k]*delta_alpha_i + Q_j[k]*delta_alpha_j; // Gk += Q[k][i]*Dai + Q[k][j]*Daj ;
	}

      // update alpha_status and G_bar

      {
	bool ui = is_upper_bound(i);
	bool uj = is_upper_bound(j);
	update_alpha_status(i);
	update_alpha_status(j);
	int k;
	if(ui != is_upper_bound(i)) // 至少有一次是upper_bound
	  {
	    Q_i = Q.get_Q(i,l);
	    if(ui) // 以前是 更新G_bar
	      for(k=0;k<l;k++)
		G_bar[k] -= C_i * Q_i[k];
	    else
	      for(k=0;k<l;k++) // 现在是 
		G_bar[k] += C_i * Q_i[k];
	  }

	if(uj != is_upper_bound(j))
	  {
	    Q_j = Q.get_Q(j,l);
	    if(uj)
	      for(k=0;k<l;k++)
		G_bar[k] -= C_j * Q_j[k];
	    else
	      for(k=0;k<l;k++)
		G_bar[k] += C_j * Q_j[k];
	  }
      }
    }

  if(iter >= max_iter)
    {
      if(active_size < l)
	{
	  // reconstruct the whole gradient to calculate objective value
	  reconstruct_gradient();
	  active_size = l;
	  info("*");
	}
      info("\nWARNING: reaching max number of iterations");
    }

  // calculate rho
  
  si->rho = calculate_rho(); // 貌似计算p松弛条件和更新b

  // calculate objective value
  // G = (Qa+p) 目标值为 1/2*a*Q*a + p*a ; 
  {
    double v = 0;
    int i;
    for(i=0;i<l;i++)
      v += alpha[i] * (G[i] + p[i]);

    si->obj = v/2;
  }

  // put back the solution
  {
    for(int i=0;i<l;i++)
      alpha_[active_set[i]] = alpha[i];
  }

  // juggle everything back
  /*{
    for(int i=0;i<l;i++)
    while(active_set[i] != i)
    swap_index(i,active_set[i]);
    // or Q.swap_index(i,active_set[i]);
    }*/

  si->upper_bound_p = Cp;
  si->upper_bound_n = Cn;

  info("\noptimization finished, #iter = %d\n",iter);

  delete[] p;
  delete[] y;
  delete[] alpha;
  delete[] alpha_status;
  delete[] active_set;
  delete[] G;
  delete[] G_bar;
}

// return 1 if already optimal, return 0 otherwise
int Solver::select_working_set(int &out_i, int &out_j)
{
  // return i,j such that
  // i: maximizes -y_i * grad(f)_i, i in I_up(\alpha)
  // j: minimizes the decrease of obj value
  //    (if quadratic coefficeint <= 0, replace it with tau)
  //    -y_j*grad(f)_j < -y_i*grad(f)_i, j in I_low(\alpha)
  /*
    KKT条件推出:
    yi = 1 , ai < C 
    -> yi*(Sum(aj*yj*kij)+b) >= 1  
    -> yi(Sum(aj*yj*kij)) + b >= 1 
    -> yi(Sum(aj*yj*kij)) -1 >= -b
    -> G(i) >= -b    b >= -G(i) ;
    同理推出：
    yi = -1 , ai > 0  ->  b >= G(i) ;
    yi = -1 , ai < C  ->  b <= G(i) ;
    yi = 1  , ai > 0  ->  b <= -G(i) ;

    通过 找寻一对 i j  使得 不满足上述不等式 
    Iup(a)  = {t| ai<C,yi= 1 or ai>0,yi=-1}
    Ilow(a) = {t| ai<C,yi=-1 or ai>0,yi= 1} 
    -yi*G(i) > -yj*G(j) when i->Iup and j->Ilow ;
  */
  double Gmax = -INF;
  double Gmax2 = -INF;
  int Gmax_idx = -1;
  int Gmin_idx = -1;
  double obj_diff_min = INF;
  // active_size 为需要来进行优化的alpha
  for(int t=0;t<active_size;t++)
    if(y[t]==+1)	
      {
	if(!is_upper_bound(t))  // yi = 1 , ai < C ;
	  if(-G[t] >= Gmax)
	    {
	      Gmax = -G[t];
	      Gmax_idx = t;
	    }
      }
    else
      {
	if(!is_lower_bound(t)) // yi = -1 , ai > 0 
	  if(G[t] >= Gmax)
	    {
	      Gmax = G[t];
	      Gmax_idx = t;
	    }
      }

  int i = Gmax_idx;
  const Qfloat *Q_i = NULL;
  if(i != -1) // NULL Q_i not accessed: Gmax=-INF if i=-1
    Q_i = Q->get_Q(i,active_size);
  // 通过计算二次梯度来继续最大化
  for(int j=0;j<active_size;j++)
    {
      if(y[j]==+1)
	{
	  if (!is_lower_bound(j)) // yj = 1 , aj > 0 
	    {
	      double grad_diff=Gmax+G[j]; // -yi*G[i]+yj*G[j] ; yi = -1 Gmax = G[i] , yi = 1 Gmax = -G[i] ; 
	      if (G[j] >= Gmax2)
		Gmax2 = G[j];
	      if (grad_diff > 0)
		{
		  double obj_diff; 
		  double quad_coef = QD[i]+QD[j]-2.0*y[i]*Q_i[j];
		  // obj_diff 为目标函数的最小值
		  if (quad_coef > 0)
		    obj_diff = -(grad_diff*grad_diff)/quad_coef;
		  else
		    obj_diff = -(grad_diff*grad_diff)/TAU;

		  if (obj_diff <= obj_diff_min)
		    {
		      Gmin_idx=j;
		      obj_diff_min = obj_diff;
		    }
		}
	    }
	}
      else
	{
	  if (!is_upper_bound(j))
	    {
	      double grad_diff= Gmax-G[j];
	      if (-G[j] >= Gmax2)
		Gmax2 = -G[j];
	      if (grad_diff > 0)
		{
		  double obj_diff; 
		  double quad_coef = QD[i]+QD[j]+2.0*y[i]*Q_i[j];
		  if (quad_coef > 0)
		    obj_diff = -(grad_diff*grad_diff)/quad_coef;
		  else
		    obj_diff = -(grad_diff*grad_diff)/TAU;

		  if (obj_diff <= obj_diff_min)
		    {
		      Gmin_idx=j;
		      obj_diff_min = obj_diff;
		    }
		}
	    }
	}
    }
  // Gmax 与 Gmax2 来统计是否迭代结束
  if(Gmax+Gmax2 < eps)
    return 1;

  out_i = Gmax_idx;
  out_j = Gmin_idx;
  return 0;
}
// 用来判断是否不用计算
bool Solver::be_shrunk(int i, double Gmax1, double Gmax2)
{
  if(is_upper_bound(i))
    {
      if(y[i]==+1)
	return(-G[i] > Gmax1);
      else
	return(-G[i] > Gmax2);
    }
  else if(is_lower_bound(i))
    {
      if(y[i]==+1)
	return(G[i] > Gmax2);
      else	
	return(G[i] > Gmax1);
    }
  else
    return(false);
}

void Solver::do_shrinking()
{
  int i;
  double Gmax1 = -INF;		// max { -y_i * grad(f)_i | i in I_up(\alpha) }
  double Gmax2 = -INF;		// max { y_i * grad(f)_i | i in I_low(\alpha) }

  // find maximal violating pair first
  for(i=0;i<active_size;i++)
    {
      if(y[i]==+1)	
	{
	  if(!is_upper_bound(i))	
	    {
	      if(-G[i] >= Gmax1)
		Gmax1 = -G[i];
	    }
	  if(!is_lower_bound(i))	
	    {
	      if(G[i] >= Gmax2)
		Gmax2 = G[i];
	    }
	}
      else	
	{
	  if(!is_upper_bound(i))	
	    {
	      if(-G[i] >= Gmax2)
		Gmax2 = -G[i];
	    }
	  if(!is_lower_bound(i))	
	    {
	      if(G[i] >= Gmax1)
		Gmax1 = G[i];
	    }
	}
    }

  if(unshrink == false && Gmax1 + Gmax2 <= eps*10) 
    {
      unshrink = true;
      reconstruct_gradient();
      active_size = l;
      info("*");
    }

  for(i=0;i<active_size;i++)
    if (be_shrunk(i, Gmax1, Gmax2))
      {
	active_size--;
	while (active_size > i)
	  {
	    if (!be_shrunk(active_size, Gmax1, Gmax2))
	      {
		swap_index(i,active_size);
		break;
	      }
	    active_size--;
	  }
      }
}

double Solver::calculate_rho()
{
  double r;
  int nr_free = 0;
  double ub = INF, lb = -INF, sum_free = 0;
  for(int i=0;i<active_size;i++)
    {
      double yG = y[i]*G[i];

      if(is_upper_bound(i))
	{
	  if(y[i]==-1)
	    ub = min(ub,yG);
	  else
	    lb = max(lb,yG);
	}
      else if(is_lower_bound(i))
	{
	  if(y[i]==+1)
	    ub = min(ub,yG);
	  else
	    lb = max(lb,yG);
	}
      else
	{
	  ++nr_free;
	  sum_free += yG;
	}
    }
  // 通过计算非支持向量的平均
  if(nr_free>0)
    r = sum_free/nr_free;
  else
    r = (ub+lb)/2;
  // else 计算支持向量中最大的两个b值 进行平均
  return r;
}
// svm 做回归
//
// Solver for nu-svm classification and regression
//
// additional constraint: e^T \alpha = constant
//
class Solver_NU : public Solver
{
public:
  Solver_NU() {}
  void Solve(int l, const QMatrix& Q, const double *p, const schar *y,
	     double *alpha, double Cp, double Cn, double eps,
	     SolutionInfo* si, int shrinking)
  {
    this->si = si;
    Solver::Solve(l,Q,p,y,alpha,Cp,Cn,eps,si,shrinking);
  }
private:
  SolutionInfo *si;
  int select_working_set(int &i, int &j);
  double calculate_rho();
  bool be_shrunk(int i, double Gmax1, double Gmax2, double Gmax3, double Gmax4);
  void do_shrinking();
};

// return 1 if already optimal, return 0 otherwise
int Solver_NU::select_working_set(int &out_i, int &out_j)
{
  // return i,j such that y_i = y_j and
  // i: maximizes -y_i * grad(f)_i, i in I_up(\alpha)
  // j: minimizes the decrease of obj value
  //    (if quadratic coefficeint <= 0, replace it with tau)
  //    -y_j*grad(f)_j < -y_i*grad(f)_i, j in I_low(\alpha)

  double Gmaxp = -INF;
  double Gmaxp2 = -INF;
  int Gmaxp_idx = -1;

  double Gmaxn = -INF;
  double Gmaxn2 = -INF;
  int Gmaxn_idx = -1;

  int Gmin_idx = -1;
  double obj_diff_min = INF;
  // 找到Iup
  for(int t=0;t<active_size;t++)
    if(y[t]==+1)
      {
	if(!is_upper_bound(t)) 
	  if(-G[t] >= Gmaxp)
	    {
	      Gmaxp = -G[t];
	      Gmaxp_idx = t;
	    }
      }
    else
      {
	if(!is_lower_bound(t))
	  if(G[t] >= Gmaxn)
	    {
	      Gmaxn = G[t];
	      Gmaxn_idx = t;
	    }
      }

  int ip = Gmaxp_idx;
  int in = Gmaxn_idx;
  const Qfloat *Q_ip = NULL;
  const Qfloat *Q_in = NULL;
  if(ip != -1) // NULL Q_ip not accessed: Gmaxp=-INF if ip=-1
    Q_ip = Q->get_Q(ip,active_size);
  if(in != -1)
    Q_in = Q->get_Q(in,active_size);

  for(int j=0;j<active_size;j++)
    {
      if(y[j]==+1)
	{
	  if (!is_lower_bound(j))	
	    {
	      double grad_diff=Gmaxp+G[j];
	      if (G[j] >= Gmaxp2)
		Gmaxp2 = G[j];
	      if (grad_diff > 0)
		{
		  double obj_diff; 
		  double quad_coef = QD[ip]+QD[j]-2*Q_ip[j];
		  if (quad_coef > 0)
		    obj_diff = -(grad_diff*grad_diff)/quad_coef;
		  else
		    obj_diff = -(grad_diff*grad_diff)/TAU;

		  if (obj_diff <= obj_diff_min)
		    {
		      Gmin_idx=j;
		      obj_diff_min = obj_diff;
		    }
		}
	    }
	}
      else
	{
	  if (!is_upper_bound(j))
	    {
	      double grad_diff=Gmaxn-G[j];
	      if (-G[j] >= Gmaxn2)
		Gmaxn2 = -G[j];
	      if (grad_diff > 0)
		{
		  double obj_diff; 
		  double quad_coef = QD[in]+QD[j]-2*Q_in[j];
		  if (quad_coef > 0)
		    obj_diff = -(grad_diff*grad_diff)/quad_coef;
		  else
		    obj_diff = -(grad_diff*grad_diff)/TAU;

		  if (obj_diff <= obj_diff_min)
		    {
		      Gmin_idx=j;
		      obj_diff_min = obj_diff;
		    }
		}
	    }
	}
    }

  if(max(Gmaxp+Gmaxp2,Gmaxn+Gmaxn2) < eps)
    return 1;

  if (y[Gmin_idx] == +1)
    out_i = Gmaxp_idx;
  else
    out_i = Gmaxn_idx;
  out_j = Gmin_idx;

  return 0;
}

bool Solver_NU::be_shrunk(int i, double Gmax1, double Gmax2, double Gmax3, double Gmax4)
{
  if(is_upper_bound(i))
    {
      if(y[i]==+1)
	return(-G[i] > Gmax1);
      else	
	return(-G[i] > Gmax4);
    }
  else if(is_lower_bound(i))
    {
      if(y[i]==+1)
	return(G[i] > Gmax2);
      else	
	return(G[i] > Gmax3);
    }
  else
    return(false);
}

void Solver_NU::do_shrinking()
{
  double Gmax1 = -INF;	// max { -y_i * grad(f)_i | y_i = +1, i in I_up(\alpha) }
  double Gmax2 = -INF;	// max { y_i * grad(f)_i | y_i = +1, i in I_low(\alpha) }
  double Gmax3 = -INF;	// max { -y_i * grad(f)_i | y_i = -1, i in I_up(\alpha) }
  double Gmax4 = -INF;	// max { y_i * grad(f)_i | y_i = -1, i in I_low(\alpha) }

  // find maximal violating pair first
  int i;
  for(i=0;i<active_size;i++)
    {
      if(!is_upper_bound(i))
	{
	  if(y[i]==+1)
	    {
	      if(-G[i] > Gmax1) Gmax1 = -G[i];
	    }
	  else	if(-G[i] > Gmax4) Gmax4 = -G[i];
	}
      if(!is_lower_bound(i))
	{
	  if(y[i]==+1)
	    {	
	      if(G[i] > Gmax2) Gmax2 = G[i];
	    }
	  else	if(G[i] > Gmax3) Gmax3 = G[i];
	}
    }

  if(unshrink == false && max(Gmax1+Gmax2,Gmax3+Gmax4) <= eps*10) 
    {
      unshrink = true;
      reconstruct_gradient();
      active_size = l;
    }

  for(i=0;i<active_size;i++)
    if (be_shrunk(i, Gmax1, Gmax2, Gmax3, Gmax4))
      {
	active_size--;
	while (active_size > i)
	  {
	    if (!be_shrunk(active_size, Gmax1, Gmax2, Gmax3, Gmax4))
	      {
		swap_index(i,active_size);
		break;
	      }
	    active_size--;
	  }
      }
}

double Solver_NU::calculate_rho()
{
  int nr_free1 = 0,nr_free2 = 0;
  double ub1 = INF, ub2 = INF;
  double lb1 = -INF, lb2 = -INF;
  double sum_free1 = 0, sum_free2 = 0;

  for(int i=0;i<active_size;i++)
    {
      if(y[i]==+1)
	{
	  if(is_upper_bound(i))
	    lb1 = max(lb1,G[i]);
	  else if(is_lower_bound(i))
	    ub1 = min(ub1,G[i]);
	  else
	    {
	      ++nr_free1;
	      sum_free1 += G[i];
	    }
	}
      else
	{
	  if(is_upper_bound(i))
	    lb2 = max(lb2,G[i]);
	  else if(is_lower_bound(i))
	    ub2 = min(ub2,G[i]);
	  else
	    {
	      ++nr_free2;
	      sum_free2 += G[i];
	    }
	}
    }

  double r1,r2;
  if(nr_free1 > 0)
    r1 = sum_free1/nr_free1;
  else
    r1 = (ub1+lb1)/2;
	
  if(nr_free2 > 0)
    r2 = sum_free2/nr_free2;
  else
    r2 = (ub2+lb2)/2;
	
  si->r = (r1+r2)/2;
  return (r1-r2)/2;
}

//
// Q matrices for various formulations
//
class SVC_Q: public Kernel
{ 
public:
  SVC_Q(const svm_problem& prob, const svm_parameter& param, const schar *y_)
    :Kernel(prob.l, prob.x, param)
  {
    clone(y,y_,prob.l);
    cache = new Cache(prob.l,(long int)(param.cache_size*(1<<20)));
    QD = new double[prob.l];
    for(int i=0;i<prob.l;i++)
      QD[i] = (this->*kernel_function)(i,i);
  }
	
  Qfloat *get_Q(int i, int len) const
  {
    Qfloat *data;
    int start, j;
    if((start = cache->get_data(i,&data,len)) < len)
      {
	for(j=start;j<len;j++)
	  data[j] = (Qfloat)(y[i]*y[j]*(this->*kernel_function)(i,j));
      }
    return data;
  }

  double *get_QD() const
  {
    return QD;
  }

  void swap_index(int i, int j) const
  {
    cache->swap_index(i,j);
    Kernel::swap_index(i,j);
    swap(y[i],y[j]);
    swap(QD[i],QD[j]);
  }

  ~SVC_Q()
  {
    delete[] y;
    delete cache;
    delete[] QD;
  }
private:
  schar *y;
  Cache *cache;
  double *QD;
};

class ONE_CLASS_Q: public Kernel
{
public:
  ONE_CLASS_Q(const svm_problem& prob, const svm_parameter& param)
    :Kernel(prob.l, prob.x, param)
  {
    cache = new Cache(prob.l,(long int)(param.cache_size*(1<<20)));
    QD = new double[prob.l];
    for(int i=0;i<prob.l;i++)
      QD[i] = (this->*kernel_function)(i,i);
  }
	
  Qfloat *get_Q(int i, int len) const
  {
    Qfloat *data;
    int start, j;
    if((start = cache->get_data(i,&data,len)) < len)
      {
	for(j=start;j<len;j++)
	  data[j] = (Qfloat)(this->*kernel_function)(i,j); // 没有标签
      }
    return data;
  }

  double *get_QD() const
  {
    return QD;
  }

  void swap_index(int i, int j) const
  {
    cache->swap_index(i,j);
    Kernel::swap_index(i,j);
    swap(QD[i],QD[j]);
  }

  ~ONE_CLASS_Q()
  {
    delete cache;
    delete[] QD;
  }
private:
  Cache *cache;
  double *QD;
};

class SVR_Q: public Kernel
{ 
public:
  SVR_Q(const svm_problem& prob, const svm_parameter& param)
    :Kernel(prob.l, prob.x, param)
  {
    l = prob.l;
    cache = new Cache(l,(long int)(param.cache_size*(1<<20)));
    QD = new double[2*l];
    sign = new schar[2*l];
    index = new int[2*l];
    for(int k=0;k<l;k++)
      {
	sign[k] = 1;
	sign[k+l] = -1;
	index[k] = k;
	index[k+l] = k;
	QD[k] = (this->*kernel_function)(k,k);
	QD[k+l] = QD[k];
      }
    buffer[0] = new Qfloat[2*l];
    buffer[1] = new Qfloat[2*l];
    next_buffer = 0;
  }

  void swap_index(int i, int j) const
  {
    swap(sign[i],sign[j]);
    swap(index[i],index[j]);
    swap(QD[i],QD[j]);
  }
	
  Qfloat *get_Q(int i, int len) const
  {
    Qfloat *data;
    int j, real_i = index[i];
    if(cache->get_data(real_i,&data,l) < l)
      {
	for(j=0;j<l;j++)
	  data[j] = (Qfloat)(this->*kernel_function)(real_i,j);
      }

    // reorder and copy
    Qfloat *buf = buffer[next_buffer];
    next_buffer = 1 - next_buffer;
    schar si = sign[i];
    for(j=0;j<len;j++)
      buf[j] = (Qfloat) si * (Qfloat) sign[j] * data[index[j]];
    return buf;
  }

  double *get_QD() const
  {
    return QD;
  }

  ~SVR_Q()
  {
    delete cache;
    delete[] sign;
    delete[] index;
    delete[] buffer[0];
    delete[] buffer[1];
    delete[] QD;
  }
private:
  int l;
  Cache *cache;
  schar *sign;
  int *index;
  mutable int next_buffer;
  Qfloat *buffer[2];
  double *QD;
};

//
// construct and solve various formulations
//
static void solve_c_svc(
			const svm_problem *prob, const svm_parameter* param,
			double *alpha, Solver::SolutionInfo* si, double Cp, double Cn)
{
  int l = prob->l;
  double *minus_ones = new double[l];
  schar *y = new schar[l];

  int i;

  for(i=0;i<l;i++)
    {
      alpha[i] = 0; // alpha = 0 保证最初条件 alpha*y成立 
      minus_ones[i] = -1; // p 数组 全为-1 
      if(prob->y[i] > 0) y[i] = +1; else y[i] = -1;
    }

  Solver s;
  s.Solve(l, SVC_Q(*prob,*param,y), minus_ones, y,
	  alpha, Cp, Cn, param->eps, si, param->shrinking); // 计算svm模型

  double sum_alpha=0;
  for(i=0;i<l;i++)
    sum_alpha += alpha[i];

  if (Cp==Cn)
    info("nu = %f\n", sum_alpha/(Cp*prob->l));

  for(i=0;i<l;i++)
    alpha[i] *= y[i];

  delete[] minus_ones;
  delete[] y;
}
//Q: what is the difference between nu-svc and c-svc?
//A: Basically they are the same thing but with different paramenters. The range of C is from zero to infinity but nu-svc is always bettween [0,1]. A nice property of nu is related to the ratio of support vectors and the ratio of the training error.
static void solve_nu_svc(const svm_problem *prob, const svm_parameter *param,
			 double *alpha, Solver::SolutionInfo* si)
{
  int i;
  int l = prob->l;
  double nu = param->nu;

  schar *y = new schar[l];

  for(i=0;i<l;i++)
    if(prob->y[i]>0)
      y[i] = +1;
    else
      y[i] = -1;

  double sum_pos = nu*l/2;
  double sum_neg = nu*l/2;

  for(i=0;i<l;i++)
    if(y[i] == +1)
      {
	alpha[i] = min(1.0,sum_pos);
	sum_pos -= alpha[i];
      }
    else
      {
	alpha[i] = min(1.0,sum_neg);
	sum_neg -= alpha[i];
      }

  double *zeros = new double[l];

  for(i=0;i<l;i++)
    zeros[i] = 0; // p is zero 

  Solver_NU s;
  s.Solve(l, SVC_Q(*prob,*param,y), zeros, y,
	  alpha, 1.0, 1.0, param->eps, si,  param->shrinking);
  double r = si->r;

  info("C = %f\n",1/r);

  for(i=0;i<l;i++)
    alpha[i] *= y[i]/r;

  si->rho /= r;
  si->obj /= (r*r);
  si->upper_bound_p = 1/r;
  si->upper_bound_n = 1/r;

  delete[] y;
  delete[] zeros;
}

static void solve_one_class(
			    const svm_problem *prob, const svm_parameter *param,
			    double *alpha, Solver::SolutionInfo* si)
{
  int l = prob->l;
  double *zeros = new double[l];
  schar *ones = new schar[l];
  int i;

  int n = (int)(param->nu*prob->l);	// # of alpha's at upper bound

  for(i=0;i<n;i++)
    alpha[i] = 1;
  if(n<prob->l)
    alpha[n] = param->nu * prob->l - n;
  for(i=n+1;i<l;i++)
    alpha[i] = 0;

  for(i=0;i<l;i++)
    {
      zeros[i] = 0;
      ones[i] = 1;
    }

  Solver s;
  s.Solve(l, ONE_CLASS_Q(*prob,*param), zeros, ones,
	  alpha, 1.0, 1.0, param->eps, si, param->shrinking);

  delete[] zeros;
  delete[] ones;
}

static void solve_epsilon_svr(
			      const svm_problem *prob, const svm_parameter *param,
			      double *alpha, Solver::SolutionInfo* si)
{
  int l = prob->l;
  double *alpha2 = new double[2*l];
  double *linear_term = new double[2*l];
  schar *y = new schar[2*l];
  int i;

  for(i=0;i<l;i++)
    {
      alpha2[i] = 0;
      linear_term[i] = param->p - prob->y[i];
      y[i] = 1;

      alpha2[i+l] = 0;
      linear_term[i+l] = param->p + prob->y[i];
      y[i+l] = -1;
    }

  Solver s;
  s.Solve(2*l, SVR_Q(*prob,*param), linear_term, y,
	  alpha2, param->C, param->C, param->eps, si, param->shrinking);

  double sum_alpha = 0;
  for(i=0;i<l;i++)
    {
      alpha[i] = alpha2[i] - alpha2[i+l];
      sum_alpha += fabs(alpha[i]);
    }
  info("nu = %f\n",sum_alpha/(param->C*l));

  delete[] alpha2;
  delete[] linear_term;
  delete[] y;
}

static void solve_nu_svr(
			 const svm_problem *prob, const svm_parameter *param,
			 double *alpha, Solver::SolutionInfo* si)
{
  int l = prob->l;
  double C = param->C;
  double *alpha2 = new double[2*l];
  double *linear_term = new double[2*l];
  schar *y = new schar[2*l];
  int i;

  double sum = C * param->nu * l / 2;
  for(i=0;i<l;i++)
    {
      alpha2[i] = alpha2[i+l] = min(sum,C);
      sum -= alpha2[i];

      linear_term[i] = - prob->y[i];
      y[i] = 1;

      linear_term[i+l] = prob->y[i];
      y[i+l] = -1;
    }

  Solver_NU s;
  s.Solve(2*l, SVR_Q(*prob,*param), linear_term, y,
	  alpha2, C, C, param->eps, si, param->shrinking);

  info("epsilon = %f\n",-si->r);

  for(i=0;i<l;i++)
    alpha[i] = alpha2[i] - alpha2[i+l];

  delete[] alpha2;
  delete[] linear_term;
  delete[] y;
}

//
// decision_function
//
struct decision_function
{
  double *alpha;
  double rho;	
};

static decision_function svm_train_one(
				       const svm_problem *prob, const svm_parameter *param,
				       double Cp, double Cn)
{
  double *alpha = Malloc(double,prob->l);
  Solver::SolutionInfo si;
  switch(param->svm_type)
    {
    case C_SVC:
      solve_c_svc(prob,param,alpha,&si,Cp,Cn);
      break;
    case NU_SVC:
      solve_nu_svc(prob,param,alpha,&si);
      break;
    case ONE_CLASS:
      solve_one_class(prob,param,alpha,&si);
      break;
    case EPSILON_SVR:
      solve_epsilon_svr(prob,param,alpha,&si);
      break;
    case NU_SVR:
      solve_nu_svr(prob,param,alpha,&si);
      break;
    }

  info("obj = %f, rho = %f\n",si.obj,si.rho);

  // output SVs

  int nSV = 0;
  int nBSV = 0;
  for(int i=0;i<prob->l;i++)
    {
      if(fabs(alpha[i]) > 0)
	{
	  ++nSV;
	  if(prob->y[i] > 0)
	    {
	      if(fabs(alpha[i]) >= si.upper_bound_p)
		++nBSV;
	    }
	  else
	    {
	      if(fabs(alpha[i]) >= si.upper_bound_n)
		++nBSV;
	    }
	}
    }

  info("nSV = %d, nBSV = %d\n",nSV,nBSV);

  decision_function f;
  f.alpha = alpha;
  f.rho = si.rho;
  return f;
}

// Platt's binary SVM Probablistic Output: an improvement from Lin et al.
// dec_values 样本的预测值
// labels 样本的目标值
// rij = 1 / (1+e^(Af+B)) 确定A B
static void sigmoid_train(
			  int l, const double *dec_values, const double *labels, 
			  double& A, double& B)
{
  double prior1=0, prior0 = 0;
  int i;

  for (i=0;i<l;i++)
    if (labels[i] > 0) prior1+=1; // 正类数统计
    else prior0+=1; // 负类数统计
	
  int max_iter=100;	// Maximal number of iterations
  double min_step=1e-10;	// Minimal step taken in line search
  double sigma=1e-12;	// For numerically strict PD of Hessian
  double eps=1e-5;
  double hiTarget=(prior1+1.0)/(prior1+2.0);// 正类数的新目标值，保证小于1，大于0
  double loTarget=1/(prior0+2.0); // 负类的新目标值
  double *t=Malloc(double,l); // 样本的新目标值数组
  double fApB,p,q,h11,h22,h21,g1,g2,det,dA,dB,gd,stepsize;
  double newA,newB,newf,d1,d2;
  int iter; 
	
  // Initial Point and Initial Fun Value
  A=0.0; B=log((prior0+1.0)/(prior1+1.0));
  double fval = 0.0;// 判断函数，越小越好，用来判断A，B的值

  for (i=0;i<l;i++)
    {
      if (labels[i]>0) t[i]=hiTarget; // 标准化判断函数的值
      else t[i]=loTarget;
      fApB = dec_values[i]*A+B;
      if (fApB>=0)
	fval += t[i]*fApB + log(1+exp(-fApB)); // 1 + exp(-fApB) = 1 / (1+e^(fApB)) ;
      else
	fval += (t[i] - 1)*fApB +log(1+exp(fApB)); // 同理
    }
  for (iter=0;iter<max_iter;iter++)
    {
      // Update Gradient and Hessian (use H' = H + sigma I)
      h11=sigma; // numerically ensures strict PD
      h22=sigma;
      h21=0.0;g1=0.0;g2=0.0;
      for (i=0;i<l;i++)
	{
	  fApB = dec_values[i]*A+B;
	  if (fApB >= 0)
	    {
	      p=exp(-fApB)/(1.0+exp(-fApB));
	      q=1.0/(1.0+exp(-fApB));
	    }
	  else
	    {
	      p=1.0/(1.0+exp(fApB));
	      q=exp(fApB)/(1.0+exp(fApB));
	    }
	  d2=p*q;
	  h11+=dec_values[i]*dec_values[i]*d2;
	  h22+=d2;
	  h21+=dec_values[i]*d2;
	  d1=t[i]-p;
	  g1+=dec_values[i]*d1;
	  g2+=d1;
	}

      // Stopping Criteria
      if (fabs(g1)<eps && fabs(g2)<eps)
	break;

      // Finding Newton direction: -inv(H') * g
      det=h11*h22-h21*h21;
      dA=-(h22*g1 - h21 * g2) / det;
      dB=-(-h21*g1+ h11 * g2) / det;
      gd=g1*dA+g2*dB;


      stepsize = 1;		// Line Search
      while (stepsize >= min_step)
	{
	  newA = A + stepsize * dA;
	  newB = B + stepsize * dB;

	  // New function value
	  newf = 0.0;
	  for (i=0;i<l;i++)
	    {
	      fApB = dec_values[i]*newA+newB;
	      if (fApB >= 0)
		newf += t[i]*fApB + log(1+exp(-fApB));
	      else
		newf += (t[i] - 1)*fApB +log(1+exp(fApB));
	    }
	  // Check sufficient decrease
	  if (newf<fval+0.0001*stepsize*gd)
	    {
	      A=newA;B=newB;fval=newf;
	      break;
	    }
	  else
	    stepsize = stepsize / 2.0;
	}

      if (stepsize < min_step)
	{
	  info("Line search fails in two-class probability estimates\n");
	  break;
	}
    }

  if (iter>=max_iter)
    info("Reaching maximal iterations in two-class probability estimates\n");
  free(t);
}
// 返回概率函数值，在函数svm_predict_probability的分类概率预测中调用，参数A，B和model的probA probB有关
static double sigmoid_predict(double decision_value, double A, double B)
{
  double fApB = decision_value*A+B;
  // 1-p used later; avoid catastrophic cancellation
  if (fApB >= 0)
    return exp(-fApB)/(1.0+exp(-fApB));
  else
    return 1.0/(1+exp(fApB)) ;
}

// Method 2 from the multiclass_prob paper by Wu, Lin, and Weng
static void multiclass_probability(int k, double **r, double *p)
{
  // k 是类的数目 
  int t,j;
  int iter = 0, max_iter=max(100,k);
  double **Q=Malloc(double *,k);
  double *Qp=Malloc(double,k);
  double pQp, eps=0.005/k;
	
  for (t=0;t<k;t++)
    {
      p[t]=1.0/k;  // Valid if k = 1
      Q[t]=Malloc(double,k);
      Q[t][t]=0;
      for (j=0;j<t;j++)
	{
	  Q[t][t]+=r[j][t]*r[j][t];
	  Q[t][j]=Q[j][t];
	}
      for (j=t+1;j<k;j++)
	{
	  Q[t][t]+=r[j][t]*r[j][t];
	  Q[t][j]=-r[j][t]*r[t][j];
	}
    }
  for (iter=0;iter<max_iter;iter++)
    {
      // stopping condition, recalculate QP,pQP for numerical accuracy
      pQp=0;
      for (t=0;t<k;t++)
	{
	  Qp[t]=0;
	  for (j=0;j<k;j++)
	    Qp[t]+=Q[t][j]*p[j];
	  pQp+=p[t]*Qp[t];
	}
      double max_error=0;
      for (t=0;t<k;t++)
	{
	  double error=fabs(Qp[t]-pQp);
	  if (error>max_error)
	    max_error=error;
	}
      if (max_error<eps) break;
		
      for (t=0;t<k;t++)
	{
	  double diff=(-Qp[t]+pQp)/Q[t][t];
	  p[t]+=diff;
	  pQp=(pQp+diff*(diff*Q[t][t]+2*Qp[t]))/(1+diff)/(1+diff);
	  for (j=0;j<k;j++)
	    {
	      Qp[j]=(Qp[j]+diff*Q[t][j])/(1+diff);
	      p[j]/=(1+diff);
	    }
	}
    }
  if (iter>=max_iter)
    info("Exceeds max_iter in multiclass_prob\n");
  for(t=0;t<k;t++) free(Q[t]);
  free(Q);
  free(Qp);
}

// Cross-validation decision values for probability estimates
/*
  k折交叉验证，将样本分为k个子样本集，然后验证其中一个样本集，用剩下的k-1个样本集为训练样本。
 */
static void svm_binary_svc_probability(
				       const svm_problem *prob, const svm_parameter *param,
				       double Cp, double Cn, double& probA, double& probB)
{
  int i;
  int nr_fold = 5;
  int *perm = Malloc(int,prob->l);
  double *dec_values = Malloc(double,prob->l);

  // random shuffle
  for(i=0;i<prob->l;i++) perm[i]=i;
  for(i=0;i<prob->l;i++)
    {
      int j = i+rand()%(prob->l-i);
      swap(perm[i],perm[j]);
    }
  for(i=0;i<nr_fold;i++)
    {
      // begin end标识测试样本集的开始与结束。
      int begin = i*prob->l/nr_fold;
      int end = (i+1)*prob->l/nr_fold;
      int j,k;
      struct svm_problem subprob;
      // 得到子训练集 k-1个子样本集。
      subprob.l = prob->l-(end-begin);
      subprob.x = Malloc(struct svm_node*,subprob.l);
      subprob.y = Malloc(double,subprob.l);
      // 填充k-1个子样本训练集
      k=0;
      for(j=0;j<begin;j++)
	{
	  subprob.x[k] = prob->x[perm[j]];
	  subprob.y[k] = prob->y[perm[j]];
	  ++k;
	}
      for(j=end;j<prob->l;j++)
	{
	  subprob.x[k] = prob->x[perm[j]];
	  subprob.y[k] = prob->y[perm[j]];
	  ++k;
	}
      int p_count=0,n_count=0;
      for(j=0;j<k;j++)
	if(subprob.y[j]>0)
	  p_count++; // 统计正样本的个数
	else
	  n_count++; // 统计负样本的个数

      if(p_count==0 && n_count==0)
	for(j=begin;j<end;j++)
	  dec_values[perm[j]] = 0;
      else if(p_count > 0 && n_count == 0)
	for(j=begin;j<end;j++)
	  dec_values[perm[j]] = 1;
      else if(p_count == 0 && n_count > 0)
	for(j=begin;j<end;j++)
	  dec_values[perm[j]] = -1;
      else
	{
	  svm_parameter subparam = *param;
	  subparam.probability=0;
	  subparam.C=1.0;
	  subparam.nr_weight=2;
	  subparam.weight_label = Malloc(int,2);
	  subparam.weight = Malloc(double,2);
	  subparam.weight_label[0]=+1;
	  subparam.weight_label[1]=-1;
	  subparam.weight[0]=Cp;
	  subparam.weight[1]=Cn;
	  struct svm_model *submodel = svm_train(&subprob,&subparam); // 训练样本集
	  for(j=begin;j<end;j++)
	    {
	      svm_predict_values(submodel,prob->x[perm[j]],&(dec_values[perm[j]])); // 预测样本集
	      // ensure +1 -1 order; reason not using CV subroutine
	      dec_values[perm[j]] *= submodel->label[0];// 标识这个预测样本是否正确  
	    }		
	  svm_free_and_destroy_model(&submodel);
	  svm_destroy_param(&subparam);
	}
      free(subprob.x);
      free(subprob.y);
    }		
  sigmoid_train(prob->l,dec_values,prob->y,probA,probB); // 计算对偶关系probA probB 
  free(dec_values);
  free(perm);
}

// Return parameter of a Laplace distribution 
// 用于回归训练 评估误差概率，返回误差参数
// 使用交叉验证
static double svm_svr_probability(
				  const svm_problem *prob, const svm_parameter *param)
{
  int i;
  int nr_fold = 5;
  double *ymv = Malloc(double,prob->l);
  double mae = 0;

  svm_parameter newparam = *param;
  newparam.probability = 0;
  svm_cross_validation(prob,&newparam,nr_fold,ymv);
  for(i=0;i<prob->l;i++)
    {
      ymv[i]=prob->y[i]-ymv[i]; // 计算交叉验证得到的标签和实际的标签的绝对值
      mae += fabs(ymv[i]);
    }		
  mae /= prob->l; // 归一化
  double std=sqrt(2*mae*mae); 
  int count=0;
  mae=0;
  for(i=0;i<prob->l;i++)
    if (fabs(ymv[i]) > 5*std) 
      count=count+1;
    else 
      mae+=fabs(ymv[i]);
  mae /= (prob->l-count);
  info("Prob. model for test data: target value = predicted value + z,\nz: Laplace distribution e^(-|z|/sigma)/(2sigma),sigma= %g\n",mae);
  free(ymv);
  return mae;
}


// label: label name, start: begin of each class, count: #data of classes, perm: indices to the original data
// perm, length l, must be allocated before calling this subroutine
// 把不同的样本目标值集中在一起
static void svm_group_classes(const svm_problem *prob, int *nr_class_ret, int **label_ret, int **start_ret, int **count_ret, int *perm)
{
  int l = prob->l;
  int max_nr_class = 16;
  int nr_class = 0;
  int *label = Malloc(int,max_nr_class); // 类目标值
  int *count = Malloc(int,max_nr_class); // 类目标值的样本数
  int *data_label = Malloc(int,l); // 数组值为该样本的类的标号 即label[data_label[i]]为该样本的真正的类目标值	
  int i;

  for(i=0;i<l;i++)
    {
      int this_label = (int)prob->y[i];
      int j;
      for(j=0;j<nr_class;j++)
	{
	  if(this_label == label[j])
	    {
	      ++count[j];
	      break;
	    }
	}
      data_label[i] = j;
      if(j == nr_class)
	{
	  if(nr_class == max_nr_class)
	    {
	      max_nr_class *= 2;
	      label = (int *)realloc(label,max_nr_class*sizeof(int));
	      count = (int *)realloc(count,max_nr_class*sizeof(int));
	    }
	  label[nr_class] = this_label;
	  count[nr_class] = 1;
	  ++nr_class;
	}
    }

  int *start = Malloc(int,nr_class); // 累加目标值的多少
  start[0] = 0;
  for(i=1;i<nr_class;i++)
    start[i] = start[i-1]+count[i-1];
  // 把不同的目标值的样本集中在不同的区段
  for(i=0;i<l;i++)
    {
      perm[start[data_label[i]]] = i;
      ++start[data_label[i]];
    }
  // 重新标识不同目标值的样本的起始位置
  start[0] = 0;
  for(i=1;i<nr_class;i++)
    start[i] = start[i-1]+count[i-1];
  
  *nr_class_ret = nr_class;
  *label_ret = label;
  *start_ret = start; // 有内存泄漏的嫌疑
  *count_ret = count;
  free(data_label);
}

//
// Interface functions
//
/*
  根据选择的算法，来组织参加训练的分类样本，以及进行训练结果的保存。其中会对样本进行初步的统计。
  1统计类别总数，同时记录类别的标签号，统计每个类的样本的数目
  2将样本按相同标签号进行分组，并连续存储
  3计算权重C
  4训练n(n-1)/2个模型，进行的是1对1的训练，以3类分类为例子
  0-1 0-2 1-2 3*（3-1）/2 个分类器 组后投票来决定多分类的标签
  5初始化nozero数组，用来统计SV(nozero sv 关联 统计支持向量的数量)
  6初始化概率数组
  7如有必要，会调用svm_binary_svc_probability
  8训练子数据集svm_train_one
  9统计一下nozero，如果nozero已经是真，就不变，如果为假，则改为真
  10 最后输出训练模型
  填充svm_model
  11 清除内存
  
  二：回归部分
  1 类别数固定为2
  2 选择性地做svm_svr_probability,one-class 不做概率评估
  3 训练
  4 输出模型
  5 清理内存
  
  分类说明：
  一对一(One-against-One,1v1)：对n类数据，每次选择其中两类，其余的类不考虑。被选择的两类数据，数据所属的类标签被重新标记为+1/-1，针对该两类训练出一个模型，这样就可以训练出n*(n-1)/2个模型。预报时，将测试样本用所有训练好的模型进行预报，，得到n(n-1)/2个预报值，根据预报值对n类进行投票，得票最多的为最终的预报值，表明测试样本的归属。
  一对多(One-Against-Rest,1VR)：对n类数据，每次选择一类，剩下的所有类作为一类。训练时，对数据所属列别标签重新标记。训练完成得到n个模型。预报时，将测试样本放入n个模型中，得到n个预报值。(如果预报值全是+1/-1的话，就麻烦了)可以根据预报值的大小来判断属于某一个类。

  训练过程函数调用：
  svm_train -> svm_train_one -> solve_c_svc
  Solver s ; 初始化s
  s.Solve(l,SVC_Q(*prob,*param,y),minus_ones,y,alpha,Cp,Cn,param->eps,si,param->shrinking);
  Solve 完成主要的工作
 */
svm_model *svm_train(const svm_problem *prob, const svm_parameter *param)
{
  svm_model *model = Malloc(svm_model,1); 
  model->param = *param;
  model->free_sv = 0;	// XXX

  if(param->svm_type == ONE_CLASS ||
     param->svm_type == EPSILON_SVR ||
     param->svm_type == NU_SVR)
    {
      // regression or one-class-svm
      model->nr_class = 2;
      model->label = NULL;
      model->nSV = NULL;
      model->probA = NULL; model->probB = NULL;
      model->sv_coef = Malloc(double *,1);

      if(param->probability && // 一类分类，不做概率估计
	 (param->svm_type == EPSILON_SVR ||
	  param->svm_type == NU_SVR))
	{// 和分类不宜样 预测时 目标值=预测值+z; z:慢住啦普拉斯分布 e^(-|z|/probA[0]) / (2probA[0]) ;
	  model->probA = Malloc(double,1);
	  model->probA[0] = svm_svr_probability(prob,param); // 5折交叉验证预测，返回啦普拉斯分布系数
	}

      decision_function f = svm_train_one(prob,param,0,0);
      model->rho = Malloc(double,1);
      model->rho[0] = f.rho;
      
      int nSV = 0;
      int i;
      // 对model进行填充
      for(i=0;i<prob->l;i++)
	if(fabs(f.alpha[i]) > 0) ++nSV;       // 统计支持向量的数量
      model->l = nSV; // model中只有支持向量
      model->SV = Malloc(svm_node *,nSV); // 支持向量的空间申请
      model->sv_coef[0] = Malloc(double,nSV);
      int j = 0;
      for(i=0;i<prob->l;i++)
	if(fabs(f.alpha[i]) > 0)
	  {
	    model->SV[j] = prob->x[i];
	    model->sv_coef[0][j] = f.alpha[i];
	    ++j;
	  }		

      free(f.alpha);
    }
  else
    {
      // classification
      int l = prob->l;
      int nr_class;
      int *label = NULL;
      int *start = NULL;
      int *count = NULL;
      int *perm = Malloc(int,l);

      // group training data of the same class
      svm_group_classes(prob,&nr_class,&label,&start,&count,perm); // 按类重新排列样本 保存在perm中。perm是一个中层数组用来指向相应的样本，实际样本的存储位置并没有改变
      if(nr_class == 1) 
	info("WARNING: training data in only one class. See README for details.\n");
		
      svm_node **x = Malloc(svm_node *,l);
      int i;
      for(i=0;i<l;i++)
	x[i] = prob->x[perm[i]];

      // calculate weighted C

      double *weighted_C = Malloc(double, nr_class);
      for(i=0;i<nr_class;i++)
	weighted_C[i] = param->C;
      for(i=0;i<param->nr_weight;i++)
	{	
	  int j;
	  for(j=0;j<nr_class;j++)
	    if(param->weight_label[i] == label[j])
	      break;
	  if(j == nr_class)
	    fprintf(stderr,"WARNING: class label %d specified in weight is not found\n", param->weight_label[i]);
	  else
	    weighted_C[j] *= param->weight[i];
	}

      // train k*(k-1)/2 models
		
      bool *nonzero = Malloc(bool,l);
      for(i=0;i<l;i++)
	nonzero[i] = false;
      decision_function *f = Malloc(decision_function,nr_class*(nr_class-1)/2);

      double *probA=NULL,*probB=NULL;
      if (param->probability)
	{
	  probA=Malloc(double,nr_class*(nr_class-1)/2);
	  probB=Malloc(double,nr_class*(nr_class-1)/2);
	}
  
      //1v1：训练多样本归属分类问题。

      int p = 0;
      for(i=0;i<nr_class;i++)
	for(int j=i+1;j<nr_class;j++)
	  {
	    svm_problem sub_prob;
	    int si = start[i], sj = start[j]; // 两个样本集的开始索引
	    int ci = count[i], cj = count[j]; // 两个样本集的样本数量
	    sub_prob.l = ci+cj; // 子问题的样本数为两个样本集的样本数亮和
	    sub_prob.x = Malloc(svm_node *,sub_prob.l);
	    sub_prob.y = Malloc(double,sub_prob.l);
	    int k;
	    // 填充子问题的样本和标签
	    for(k=0;k<ci;k++)
	      {
		sub_prob.x[k] = x[si+k];
		sub_prob.y[k] = +1;
	      }
	    for(k=0;k<cj;k++)
	      {
		sub_prob.x[ci+k] = x[sj+k];
		sub_prob.y[ci+k] = -1;
	      }
	    // 预测对偶概率 得到probA probB
	    if(param->probability)
	      svm_binary_svc_probability(&sub_prob,param,weighted_C[i],weighted_C[j],probA[p],probB[p]);

	    f[p] = svm_train_one(&sub_prob,param,weighted_C[i],weighted_C[j]); // 当前模型的训练 并得到alpha和b
	    for(k=0;k<ci;k++)
	      if(!nonzero[si+k] && fabs(f[p].alpha[k]) > 0) // 统计支持向量
		nonzero[si+k] = true;
	    for(k=0;k<cj;k++)
	      if(!nonzero[sj+k] && fabs(f[p].alpha[ci+k]) > 0) 
		nonzero[sj+k] = true;
	    free(sub_prob.x);
	    free(sub_prob.y);
	    ++p;
	  }

      // build output 输出model

      model->nr_class = nr_class;
		
      model->label = Malloc(int,nr_class);
      for(i=0;i<nr_class;i++)
	model->label[i] = label[i];
		
      model->rho = Malloc(double,nr_class*(nr_class-1)/2);
      for(i=0;i<nr_class*(nr_class-1)/2;i++)
	model->rho[i] = f[i].rho;

      if(param->probability) // 保存对偶概率。
	{
	  model->probA = Malloc(double,nr_class*(nr_class-1)/2);
	  model->probB = Malloc(double,nr_class*(nr_class-1)/2);
	  for(i=0;i<nr_class*(nr_class-1)/2;i++)
	    {
	      model->probA[i] = probA[i];
	      model->probB[i] = probB[i];
	    }
	}
      else
	{
	  model->probA=NULL;
	  model->probB=NULL;
	}

      int total_sv = 0;
      int *nz_count = Malloc(int,nr_class);
      model->nSV = Malloc(int,nr_class);
      // 统计所有的支持向量
      for(i=0;i<nr_class;i++)
	{
	  int nSV = 0;
	  for(int j=0;j<count[i];j++)
	    if(nonzero[start[i]+j])
	      {	
		++nSV;
		++total_sv;
	      }
	  model->nSV[i] = nSV;
	  nz_count[i] = nSV;
	}
		
      info("Total nSV = %d\n",total_sv);

      model->l = total_sv;
      model->SV = Malloc(svm_node *,total_sv);
      p = 0;
      for(i=0;i<l;i++)
	if(nonzero[i]) model->SV[p++] = x[i];

      int *nz_start = Malloc(int,nr_class);
      nz_start[0] = 0;
      // 统计各个类别的支持向量的数量 然后进行累加
      for(i=1;i<nr_class;i++)
	nz_start[i] = nz_start[i-1]+nz_count[i-1];

      model->sv_coef = Malloc(double *,nr_class-1);
      for(i=0;i<nr_class-1;i++)
	model->sv_coef[i] = Malloc(double,total_sv);

      p = 0;
      for(i=0;i<nr_class;i++)
	for(int j=i+1;j<nr_class;j++)
	  {
	    // classifier (i,j): coefficients with
	    // i are in sv_coef[j-1][nz_start[i]...],
	    // j are in sv_coef[i][nz_start[j]...]
	    // 排列?????
	    int si = start[i];
	    int sj = start[j];
	    int ci = count[i];
	    int cj = count[j];
				
	    int q = nz_start[i];
	    int k;
	    for(k=0;k<ci;k++)
	      if(nonzero[si+k])
		model->sv_coef[j-1][q++] = f[p].alpha[k];
	    q = nz_start[j];
	    for(k=0;k<cj;k++)
	      if(nonzero[sj+k])
		model->sv_coef[i][q++] = f[p].alpha[ci+k];
	    ++p;
	  }
		
      free(label);
      free(probA);
      free(probB);
      free(count);
      free(perm);
      free(start);
      free(x);
      free(weighted_C);
      free(nonzero);
      for(i=0;i<nr_class*(nr_class-1)/2;i++)
	free(f[i].alpha);
      free(f);
      free(nz_count);
      free(nz_start);
    }
  return model;
}

// Stratified cross validation
/*
  先随机打乱次序，然后根据n折的数目，留一份作为测试集，其他的作为训练集，做n次。
 */
void svm_cross_validation(const svm_problem *prob, const svm_parameter *param, int nr_fold, double *target)
{
  int i;
  int *fold_start = Malloc(int,nr_fold+1);
  int l = prob->l;
  int *perm = Malloc(int,l);
  int nr_class;

  // stratified cv may not give leave-one-out rate
  // Each class to l folds -> some folds may have zero elements
  if((param->svm_type == C_SVC ||
      param->svm_type == NU_SVC) && nr_fold < l)
    {
      int *start = NULL;
      int *label = NULL;
      int *count = NULL;
      svm_group_classes(prob,&nr_class,&label,&start,&count,perm);// 重新排列多样本到perm中

      // random shuffle and then data grouped by fold using the array perm
      int *fold_count = Malloc(int,nr_fold);
      int c;
      int *index = Malloc(int,l);
      for(i=0;i<l;i++)
	index[i]=perm[i];
      for (c=0; c<nr_class; c++) 
	for(i=0;i<count[c];i++)
	  {
	    int j = i+rand()%(count[c]-i);
	    swap(index[start[c]+j],index[start[c]+i]);
	  }
      for(i=0;i<nr_fold;i++)
	{
	  fold_count[i] = 0;
	  for (c=0; c<nr_class;c++)
	    fold_count[i]+=(i+1)*count[c]/nr_fold-i*count[c]/nr_fold; 
	}
      fold_start[0]=0;
      for (i=1;i<=nr_fold;i++)
	fold_start[i] = fold_start[i-1]+fold_count[i-1];
      for (c=0; c<nr_class;c++)
	for(i=0;i<nr_fold;i++)
	  {
	    int begin = start[c]+i*count[c]/nr_fold;
	    int end = start[c]+(i+1)*count[c]/nr_fold;
	    for(int j=begin;j<end;j++)
	      {
		perm[fold_start[i]] = index[j];
		fold_start[i]++;
	      }
	  }
      fold_start[0]=0;
      for (i=1;i<=nr_fold;i++)
	fold_start[i] = fold_start[i-1]+fold_count[i-1];
      free(start);	
      free(label);
      free(count);	
      free(index);
      free(fold_count);
    }
  else // 回归
    {
      for(i=0;i<l;i++) perm[i]=i;
      for(i=0;i<l;i++)
	{
	  int j = i+rand()%(l-i);
	  swap(perm[i],perm[j]);
	}
      for(i=0;i<=nr_fold;i++)
	fold_start[i]=i*l/nr_fold;
    }

  for(i=0;i<nr_fold;i++)
    {
      int begin = fold_start[i];
      int end = fold_start[i+1];
      int j,k;
      struct svm_problem subprob;

      subprob.l = l-(end-begin);
      subprob.x = Malloc(struct svm_node*,subprob.l);
      subprob.y = Malloc(double,subprob.l);
			
      k=0;
      for(j=0;j<begin;j++)
	{
	  subprob.x[k] = prob->x[perm[j]];
	  subprob.y[k] = prob->y[perm[j]];
	  ++k;
	}
      for(j=end;j<l;j++)
	{
	  subprob.x[k] = prob->x[perm[j]];
	  subprob.y[k] = prob->y[perm[j]];
	  ++k;
	}
      struct svm_model *submodel = svm_train(&subprob,param);
      if(param->probability && 
	 (param->svm_type == C_SVC || param->svm_type == NU_SVC))
	{
	  double *prob_estimates=Malloc(double,svm_get_nr_class(submodel));
	  for(j=begin;j<end;j++)
	    target[perm[j]] = svm_predict_probability(submodel,prob->x[perm[j]],prob_estimates);
	  free(prob_estimates);			
	}
      else
	for(j=begin;j<end;j++)
	  target[perm[j]] = svm_predict(submodel,prob->x[perm[j]]);
      svm_free_and_destroy_model(&submodel);
      free(subprob.x);
      free(subprob.y);
    }		
  free(fold_start);
  free(perm);	
}

// get svm_type
int svm_get_svm_type(const svm_model *model)
{
  return model->param.svm_type;
}
// 得到样本的种类数
int svm_get_nr_class(const svm_model *model)
{
  return model->nr_class;
}

void svm_get_labels(const svm_model *model, int* label)
{
  if (model->label != NULL)
    for(int i=0;i<model->nr_class;i++)
      label[i] = model->label[i];
}

double svm_get_svr_probability(const svm_model *model)
{// 返回训练好的模型的概率
  if ((model->param.svm_type == EPSILON_SVR || model->param.svm_type == NU_SVR) &&
      model->probA!=NULL)
    return model->probA[0];
  else
    {
      fprintf(stderr,"Model doesn't contain information for SVR probability inference\n");
      return 0;
    }
}
// 预测样本数据的目标值，保留在double * dec_values中，
// 如果是做分类的问题，返回一堆值，供后续函数做决策；如果是回归就返回一个值。
// 其中one-v-one方法需要n(n-1)/2次，产生n(n-1)/2个预报值。
double svm_predict_values(const svm_model *model, const svm_node *x, double* dec_values)
{
  int i;
  if(model->param.svm_type == ONE_CLASS ||
     model->param.svm_type == EPSILON_SVR ||
     model->param.svm_type == NU_SVR)
    {
      // 分类问题，单个聚类
      double *sv_coef = model->sv_coef[0];
      double sum = 0;
      for(i=0;i<model->l;i++)
	sum += sv_coef[i] * Kernel::k_function(x,model->SV[i],model->param);
      sum -= model->rho[0];
      *dec_values = sum;

      if(model->param.svm_type == ONE_CLASS)
	return (sum>0)?1:-1;
      else
	return sum;
    }
  else
    {
      int nr_class = model->nr_class;
      int l = model->l;
		
      double *kvalue = Malloc(double,l);
      for(i=0;i<l;i++)
	kvalue[i] = Kernel::k_function(x,model->SV[i],model->param);

      int *start = Malloc(int,nr_class);
      start[0] = 0;
      for(i=1;i<nr_class;i++)
	start[i] = start[i-1]+model->nSV[i-1];

      int *vote = Malloc(int,nr_class);
      for(i=0;i<nr_class;i++)
	vote[i] = 0;

      int p=0;
      for(i=0;i<nr_class;i++)
	for(int j=i+1;j<nr_class;j++)
	  {
	    double sum = 0;
	    int si = start[i];
	    int sj = start[j];
	    int ci = model->nSV[i];
	    int cj = model->nSV[j];
				
	    int k;
	    double *coef1 = model->sv_coef[j-1];
	    double *coef2 = model->sv_coef[i];
	    for(k=0;k<ci;k++)
	      sum += coef1[si+k] * kvalue[si+k];
	    for(k=0;k<cj;k++)
	      sum += coef2[sj+k] * kvalue[sj+k];
	    sum -= model->rho[p];
	    dec_values[p] = sum;

	    if(dec_values[p] > 0)
	      ++vote[i];
	    else
	      ++vote[j];
	    p++;
	  }

      int vote_max_idx = 0;
      for(i=1;i<nr_class;i++)
	if(vote[i] > vote[vote_max_idx])
	  vote_max_idx = i;

      free(kvalue);
      free(start);
      free(vote);
      return model->label[vote_max_idx];
    }
}

double svm_predict(const svm_model *model, const svm_node *x)
{
  int nr_class = model->nr_class;
  double *dec_values;
  if(model->param.svm_type == ONE_CLASS ||
     model->param.svm_type == EPSILON_SVR ||
     model->param.svm_type == NU_SVR)
    dec_values = Malloc(double, 1);
  else 
    dec_values = Malloc(double, nr_class*(nr_class-1)/2);
  double pred_result = svm_predict_values(model, x, dec_values);
  free(dec_values);
  return pred_result;
}
/*
  根据model文件中对偶概率和预测值 预测各类的概率 返回最大概率的类 也就是用概率预测样本的目标值
 */
double svm_predict_probability(
			       const svm_model *model, const svm_node *x, double *prob_estimates)
{
  if ((model->param.svm_type == C_SVC || model->param.svm_type == NU_SVC) &&
      model->probA!=NULL && model->probB!=NULL)
    {
      int i;
      int nr_class = model->nr_class;
      double *dec_values = Malloc(double, nr_class*(nr_class-1)/2);
      svm_predict_values(model, x, dec_values);

      double min_prob=1e-7;
      double **pairwise_prob=Malloc(double *,nr_class);
      for(i=0;i<nr_class;i++)
	pairwise_prob[i]=Malloc(double,nr_class);
      int k=0;
      for(i=0;i<nr_class;i++)
	for(int j=i+1;j<nr_class;j++)
	  {
	    pairwise_prob[i][j]=min(max(sigmoid_predict(dec_values[k],model->probA[k],model->probB[k]),min_prob),1-min_prob);
	    pairwise_prob[j][i]=1-pairwise_prob[i][j];
	    k++;
	  }
      multiclass_probability(nr_class,pairwise_prob,prob_estimates); // prob_estimate 返回各类的概率数值，初始为0

      int prob_max_idx = 0;
      for(i=1;i<nr_class;i++)
	if(prob_estimates[i] > prob_estimates[prob_max_idx])
	  prob_max_idx = i;
      for(i=0;i<nr_class;i++)
	free(pairwise_prob[i]);
      free(dec_values);
      free(pairwise_prob);	     
      return model->label[prob_max_idx];
    }
  else 
    return svm_predict(model, x);
}

static const char *svm_type_table[] =
  {
    "c_svc","nu_svc","one_class","epsilon_svr","nu_svr",NULL
  };

static const char *kernel_type_table[]=
  {
    "linear","polynomial","rbf","sigmoid","precomputed",NULL
  };

int svm_save_model(const char *model_file_name, const svm_model *model)
{
  FILE *fp = fopen(model_file_name,"w");
  if(fp==NULL) return -1;

  char *old_locale = strdup(setlocale(LC_ALL, NULL));
  setlocale(LC_ALL, "C");

  const svm_parameter& param = model->param;

  fprintf(fp,"svm_type %s\n", svm_type_table[param.svm_type]);
  fprintf(fp,"kernel_type %s\n", kernel_type_table[param.kernel_type]);

  if(param.kernel_type == POLY)
    fprintf(fp,"degree %d\n", param.degree);

  if(param.kernel_type == POLY || param.kernel_type == RBF || param.kernel_type == SIGMOID)
    fprintf(fp,"gamma %g\n", param.gamma);

  if(param.kernel_type == POLY || param.kernel_type == SIGMOID)
    fprintf(fp,"coef0 %g\n", param.coef0);

  int nr_class = model->nr_class;
  int l = model->l;
  fprintf(fp, "nr_class %d\n", nr_class);
  fprintf(fp, "total_sv %d\n",l);
	
  {
    fprintf(fp, "rho");
    for(int i=0;i<nr_class*(nr_class-1)/2;i++)
      fprintf(fp," %g",model->rho[i]);
    fprintf(fp, "\n");
  }
	
  if(model->label)
    {
      fprintf(fp, "label");
      for(int i=0;i<nr_class;i++)
	fprintf(fp," %d",model->label[i]);
      fprintf(fp, "\n");
    }

  if(model->probA) // regression has probA only
    {
      fprintf(fp, "probA");
      for(int i=0;i<nr_class*(nr_class-1)/2;i++)
	fprintf(fp," %g",model->probA[i]);
      fprintf(fp, "\n");
    }
  if(model->probB)
    {
      fprintf(fp, "probB");
      for(int i=0;i<nr_class*(nr_class-1)/2;i++)
	fprintf(fp," %g",model->probB[i]);
      fprintf(fp, "\n");
    }

  if(model->nSV)
    {
      fprintf(fp, "nr_sv");
      for(int i=0;i<nr_class;i++)
	fprintf(fp," %d",model->nSV[i]);
      fprintf(fp, "\n");
    }

  fprintf(fp, "SV\n");
  const double * const *sv_coef = model->sv_coef;
  const svm_node * const *SV = model->SV;

  for(int i=0;i<l;i++)
    {
      for(int j=0;j<nr_class-1;j++)
	fprintf(fp, "%.16g ",sv_coef[j][i]);

      const svm_node *p = SV[i];

      if(param.kernel_type == PRECOMPUTED)
	fprintf(fp,"0:%d ",(int)(p->value));
      else
	while(p->index != -1)
	  {
	    fprintf(fp,"%d:%.8g ",p->index,p->value);
	    p++;
	  }
      fprintf(fp, "\n");
    }

  setlocale(LC_ALL, old_locale);
  free(old_locale);

  if (ferror(fp) != 0 || fclose(fp) != 0) return -1;
  else return 0;
}

static char *line = NULL;
static int max_line_len;

static char* readline(FILE *input)
{
  int len;

  if(fgets(line,max_line_len,input) == NULL)
    return NULL;

  while(strrchr(line,'\n') == NULL)
    {
      max_line_len *= 2;
      line = (char *) realloc(line,max_line_len);
      len = (int) strlen(line);
      if(fgets(line+len,max_line_len-len,input) == NULL)
	break;
    }
  return line;
}

svm_model *svm_load_model(const char *model_file_name)
{
  FILE *fp = fopen(model_file_name,"rb");
  if(fp==NULL) return NULL;

  char *old_locale = strdup(setlocale(LC_ALL, NULL));
  setlocale(LC_ALL, "C");

  // read parameters

  svm_model *model = Malloc(svm_model,1);
  svm_parameter& param = model->param;
  model->rho = NULL;
  model->probA = NULL;
  model->probB = NULL;
  model->label = NULL;
  model->nSV = NULL;

  char cmd[81];
  while(1)
    {
      fscanf(fp,"%80s",cmd);

      if(strcmp(cmd,"svm_type")==0)
	{
	  fscanf(fp,"%80s",cmd);
	  int i;
	  for(i=0;svm_type_table[i];i++)
	    {
	      if(strcmp(svm_type_table[i],cmd)==0)
		{
		  param.svm_type=i;
		  break;
		}
	    }
	  if(svm_type_table[i] == NULL)
	    {
	      fprintf(stderr,"unknown svm type.\n");
				
	      setlocale(LC_ALL, old_locale);
	      free(old_locale);
	      free(model->rho);
	      free(model->label);
	      free(model->nSV);
	      free(model);
	      return NULL;
	    }
	}
      else if(strcmp(cmd,"kernel_type")==0)
	{		
	  fscanf(fp,"%80s",cmd);
	  int i;
	  for(i=0;kernel_type_table[i];i++)
	    {
	      if(strcmp(kernel_type_table[i],cmd)==0)
		{
		  param.kernel_type=i;
		  break;
		}
	    }
	  if(kernel_type_table[i] == NULL)
	    {
	      fprintf(stderr,"unknown kernel function.\n");
				
	      setlocale(LC_ALL, old_locale);
	      free(old_locale);
	      free(model->rho);
	      free(model->label);
	      free(model->nSV);
	      free(model);
	      return NULL;
	    }
	}
      else if(strcmp(cmd,"degree")==0)
	fscanf(fp,"%d",&param.degree);
      else if(strcmp(cmd,"gamma")==0)
	fscanf(fp,"%lf",&param.gamma);
      else if(strcmp(cmd,"coef0")==0)
	fscanf(fp,"%lf",&param.coef0);
      else if(strcmp(cmd,"nr_class")==0)
	fscanf(fp,"%d",&model->nr_class);
      else if(strcmp(cmd,"total_sv")==0)
	fscanf(fp,"%d",&model->l);
      else if(strcmp(cmd,"rho")==0)
	{
	  int n = model->nr_class * (model->nr_class-1)/2;
	  model->rho = Malloc(double,n);
	  for(int i=0;i<n;i++)
	    fscanf(fp,"%lf",&model->rho[i]);
	}
      else if(strcmp(cmd,"label")==0)
	{
	  int n = model->nr_class;
	  model->label = Malloc(int,n);
	  for(int i=0;i<n;i++)
	    fscanf(fp,"%d",&model->label[i]);
	}
      else if(strcmp(cmd,"probA")==0)
	{
	  int n = model->nr_class * (model->nr_class-1)/2;
	  model->probA = Malloc(double,n);
	  for(int i=0;i<n;i++)
	    fscanf(fp,"%lf",&model->probA[i]);
	}
      else if(strcmp(cmd,"probB")==0)
	{
	  int n = model->nr_class * (model->nr_class-1)/2;
	  model->probB = Malloc(double,n);
	  for(int i=0;i<n;i++)
	    fscanf(fp,"%lf",&model->probB[i]);
	}
      else if(strcmp(cmd,"nr_sv")==0)
	{
	  int n = model->nr_class;
	  model->nSV = Malloc(int,n);
	  for(int i=0;i<n;i++)
	    fscanf(fp,"%d",&model->nSV[i]);
	}
      else if(strcmp(cmd,"SV")==0)
	{
	  while(1)
	    {
	      int c = getc(fp);
	      if(c==EOF || c=='\n') break;	
	    }
	  break;
	}
      else
	{
	  fprintf(stderr,"unknown text in model file: [%s]\n",cmd);
			
	  setlocale(LC_ALL, old_locale);
	  free(old_locale);
	  free(model->rho);
	  free(model->label);
	  free(model->nSV);
	  free(model);
	  return NULL;
	}
    }

  // read sv_coef and SV

  int elements = 0;
  long pos = ftell(fp);

  max_line_len = 1024;
  line = Malloc(char,max_line_len);
  char *p,*endptr,*idx,*val;

  while(readline(fp)!=NULL)
    {
      p = strtok(line,":");
      while(1)
	{
	  p = strtok(NULL,":");
	  if(p == NULL)
	    break;
	  ++elements;
	}
    }
  elements += model->l;

  fseek(fp,pos,SEEK_SET);

  int m = model->nr_class - 1;
  int l = model->l;
  model->sv_coef = Malloc(double *,m);
  int i;
  for(i=0;i<m;i++)
    model->sv_coef[i] = Malloc(double,l);
  model->SV = Malloc(svm_node*,l);
  svm_node *x_space = NULL;
  if(l>0) x_space = Malloc(svm_node,elements);

  int j=0;
  for(i=0;i<l;i++)
    {
      readline(fp);
      model->SV[i] = &x_space[j];

      p = strtok(line, " \t");
      model->sv_coef[0][i] = strtod(p,&endptr);
      for(int k=1;k<m;k++)
	{
	  p = strtok(NULL, " \t");
	  model->sv_coef[k][i] = strtod(p,&endptr);
	}

      while(1)
	{
	  idx = strtok(NULL, ":");
	  val = strtok(NULL, " \t");

	  if(val == NULL)
	    break;
	  x_space[j].index = (int) strtol(idx,&endptr,10);
	  x_space[j].value = strtod(val,&endptr);

	  ++j;
	}
      x_space[j++].index = -1;
    }
  free(line);

  setlocale(LC_ALL, old_locale);
  free(old_locale);

  if (ferror(fp) != 0 || fclose(fp) != 0)
    return NULL;

  model->free_sv = 1;	// XXX
  return model;
}

void svm_free_model_content(svm_model* model_ptr)
{
  if(model_ptr->free_sv && model_ptr->l > 0 && model_ptr->SV != NULL)
    free((void *)(model_ptr->SV[0]));
  if(model_ptr->sv_coef)
    {
      for(int i=0;i<model_ptr->nr_class-1;i++)
	free(model_ptr->sv_coef[i]);
    }

  free(model_ptr->SV);
  model_ptr->SV = NULL;

  free(model_ptr->sv_coef);
  model_ptr->sv_coef = NULL;

  free(model_ptr->rho);
  model_ptr->rho = NULL;

  free(model_ptr->label);
  model_ptr->label= NULL;

  free(model_ptr->probA);
  model_ptr->probA = NULL;

  free(model_ptr->probB);
  model_ptr->probB= NULL;

  free(model_ptr->nSV);
  model_ptr->nSV = NULL;
}

void svm_free_and_destroy_model(svm_model** model_ptr_ptr)
{
  if(model_ptr_ptr != NULL && *model_ptr_ptr != NULL)
    {
      svm_free_model_content(*model_ptr_ptr);
      free(*model_ptr_ptr);
      *model_ptr_ptr = NULL;
    }
}

void svm_destroy_param(svm_parameter* param)
{
  free(param->weight_label);
  free(param->weight);
}

const char *svm_check_parameter(const svm_problem *prob, const svm_parameter *param)
{
  // svm_type

  int svm_type = param->svm_type;
  if(svm_type != C_SVC &&
     svm_type != NU_SVC &&
     svm_type != ONE_CLASS &&
     svm_type != EPSILON_SVR &&
     svm_type != NU_SVR)
    return "unknown svm type";
	
  // kernel_type, degree
	
  int kernel_type = param->kernel_type;
  if(kernel_type != LINEAR &&
     kernel_type != POLY &&
     kernel_type != RBF &&
     kernel_type != SIGMOID &&
     kernel_type != PRECOMPUTED)
    return "unknown kernel type";

  if(param->gamma < 0)
    return "gamma < 0";

  if(param->degree < 0)
    return "degree of polynomial kernel < 0";

  // cache_size,eps,C,nu,p,shrinking

  if(param->cache_size <= 0)
    return "cache_size <= 0";

  if(param->eps <= 0)
    return "eps <= 0";

  if(svm_type == C_SVC ||
     svm_type == EPSILON_SVR ||
     svm_type == NU_SVR)
    if(param->C <= 0)
      return "C <= 0";

  if(svm_type == NU_SVC ||
     svm_type == ONE_CLASS ||
     svm_type == NU_SVR)
    if(param->nu <= 0 || param->nu > 1)
      return "nu <= 0 or nu > 1";

  if(svm_type == EPSILON_SVR)
    if(param->p < 0)
      return "p < 0";

  if(param->shrinking != 0 &&
     param->shrinking != 1)
    return "shrinking != 0 and shrinking != 1";

  if(param->probability != 0 &&
     param->probability != 1)
    return "probability != 0 and probability != 1";

  if(param->probability == 1 &&
     svm_type == ONE_CLASS)
    return "one-class SVM probability output not supported yet";


  // check whether nu-svc is feasible
	
  if(svm_type == NU_SVC)
    {
      int l = prob->l;
      int max_nr_class = 16;
      int nr_class = 0;
      int *label = Malloc(int,max_nr_class);
      int *count = Malloc(int,max_nr_class);

      int i;
      for(i=0;i<l;i++)
	{
	  int this_label = (int)prob->y[i];
	  int j;
	  for(j=0;j<nr_class;j++)
	    if(this_label == label[j])
	      {
		++count[j];
		break;
	      }
	  if(j == nr_class)
	    {
	      if(nr_class == max_nr_class)
		{
		  max_nr_class *= 2;
		  label = (int *)realloc(label,max_nr_class*sizeof(int));
		  count = (int *)realloc(count,max_nr_class*sizeof(int));
		}
	      label[nr_class] = this_label;
	      count[nr_class] = 1;
	      ++nr_class;
	    }
	}
	
      for(i=0;i<nr_class;i++)
	{
	  int n1 = count[i];
	  for(int j=i+1;j<nr_class;j++)
	    {
	      int n2 = count[j];
	      if(param->nu*(n1+n2)/2 > min(n1,n2))
		{
		  free(label);
		  free(count);
		  return "specified nu is infeasible";
		}
	    }
	}
      free(label);
      free(count);
    }

  return NULL;
}

int svm_check_probability_model(const svm_model *model)
{
  return ((model->param.svm_type == C_SVC || model->param.svm_type == NU_SVC) &&
	  model->probA!=NULL && model->probB!=NULL) ||
    ((model->param.svm_type == EPSILON_SVR || model->param.svm_type == NU_SVR) &&
     model->probA!=NULL);
}

void svm_set_print_string_function(void (*print_func)(const char *))
{
  if(print_func == NULL)
    svm_print_string = &print_string_stdout;
  else
    svm_print_string = print_func;
}
